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Abstract 

We propose and study an improved method to calculate the fermionic determinant of dynamical 



^ . configurations. The evaluation or at least stochastic estimation of ratios of fermionic determinants 

(N . ^ 

is essential for a recently proposed updating method of smeared link dynamical fermions. This 

^ ■ update creates a sequence of configurations by changing a subset of the gauge links by a pure gauge 

Q . heat bath or over relaxation step. The acceptance of the proposed configuration depends on the 

cn , . 

("^ , ratio of the fermionic determinants on the new and original configurations. We study this ratio as 
(N 

O ' the function of the number of links that are changed in the heat bath update. We find that even 

■ when every link of a given direction and parity of a lOfm configuration is updated, the average of 

Qh, the determinant ratio is still close to one and with the improved stochastic estimator the proposed 



change is accepted with about 20% probability. This improvement suggests that the new updating 
technique can be efficient even on large lattices. 
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I. INTRODUCTION 



The use of smeared or fat links in staggered fermion actions has gained popularity in 
recent years due to the improved flavor symmetry these actions possess ||, |[. Smeared 
links are naturally part of most systematic improvement programs and many overlap fermion 
formulations as well H, |^. The main difficulty that limits the use of smeared link fermions 
is their potential complexity in dynamical simulations. Unless the smeared links are lin- 
ear combinations of the original thin links the explicit form of the fermionic force needed 
for standard molecular dynamics simulations is very complicated, making the HMC or R 
algorithms impractical or even impossible. A recently proposed update for smeared link 
dynamical fermions |||, |^ avoids this problem by creating a sequence of configurations by 
updating a subset of the gauge links by a pure gauge heat bath or over relaxation step. 
The proposed configuration is accepted or rejected according to the change in the fermionic 
determinant. In fact one does not even have to evaluate the change in the determinant, a 
stochastic estimator can be used instead. That requires no more than the evaluation of the 
inverse fermion matrix on a Gaussian random source vector. 

The above outlined algorithm can fail in two ways. First, if the ratio of the fermionic de- 
terminants is small, the acceptance rate is small. The algorithm can also fail if the stochastic 
estimator gives a poor approximation of the determinant, making the autocorrelation time 
of the simulation (i.e. the number of independent Gaussian random sources needed to get 
a good estimate for the determinant) very large. In this paper we discuss systematic ways 
to improve the stochastic estimator. With the improved estimator we calculate the ratio 
of fermionic determinants as the function of the number of links updated with a heat bath 
step, and show that it remains close to one even if the updated volume is large. We illustrate 
and test the method using staggered fermions with HYP smeared gauge links though the 
generalization to any other smeared link action is straightforward. 

II. THE HYP ACTION AND ITS DYNAMICAL UPDATE 

In this section we define the HYP action and briefly summarize the partial-global updating 
technique. We consider a smeared link action of the form 

S = S,{U) + S,{V) + Sf{V) (1) 
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where Sg{U) and Sg(y) are gauge actions depending on the thin links {U} and smeared 
links {V}, respectively, and Sf is the fermionic action depending on the smeared links only. 
The updating method and all its improvements that we discuss in this paper would work 
with any kind of smeared links {V}, though the efficiency suffers if the smeared links are 
not smooth enough. In our work we use HYP smeared links with staggered fermions. The 
HYP links are optimized non-perturbatively to be maximally smooth. The construction and 
properties of HYP smearing are discussed in detail in Ref. |[l|]. 
We use a plaquette gauge action for Sg{U) 

5,(f/) = -f EReTr(f/,). (2) 

We choose SgiV) to improve computational efficiency and we will discuss our specific choice 
in Sect. III.C. The staggered fermionic matrix is defined in the usual way 

The matrix M\V)M(y) is block diagonal on even and odd lattice sites. In the following 
we will denote the even block by Q 

QiV) = {M\V)M{V))e.en,even (4) 

and define the fermionic action as 

Sf{V) = -^TT\nn{V) (5) 

to describe nj flavors of staggered fermions. In the following we consider Uf = 4 flavors but 
we will briefly describe the generalization to arbitrary flavors at the end of Sect. III.B. 

In Refs. [|[ 1^ a partial-global heat bath and over relaxation updating method was 
proposed to simulate the system described by Eq. |I|. In this paper we are not concerned 
about the update itself, but to motivate our interest in calculating the fermionic determinant 
ratios we briefly summarize the main points of the method. In the first step of the update 
one changes a subset of the thin links {U} to propose a new thin gauge link configuration 
{[/'}. The new links are chosen with a heat bath or over relaxed update that satisfies the 
detailed balance condition with the thin link gauge action Sg{U). The smeared links {V^} 
and {V^'} are unique once the thin links are defined. Next the proposed configuration is 
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accepted with the probability 

P^^ = min{l, exp(-S,(V') + Sg(V))^^^}. (6) 

The ratio of the determinants can be written as 

dct{Q{V')) ^ J dgeexp(-rn-i(yon(y)g) 

det{Q{V)) /<^seexp(-e«) 

= < eM-C[^-'iv')n{v) - 1]0 • (7) 

The expectation value can be evaluated stochastically where on every gauge configuration 
pair {U} and {U'} only one random source ^ is used to estimate the determinant ratio and 
the expectation value is taken together with the configuration ensemble average. That leads 
to the stochastic acceptance probability 



III. IMPROVING THE PARTIAL-GLOBAL UPDATE 

The success of the partial-global updating algorithm depends on two things. First, on 
the ratio of the determinants of the new and old links, and next on the effectiveness of 
the stochastic estimator. If the stochastic estimator fluctuates wildly, it can reduce the 
acceptance rate to practically zero even if the change in the determinant is actually small. 
In the following we will discuss improving the stochastic estimator first. 



A. Improving the stochastic estimator 

To calculate the acceptance probability we have to calculate the ratio of the determinants 

det fl' 

det -\A) = =< exp(-r[^ - 1]0 >e?=< exp(-A5/) > (9) 

where fl — fl{V), fl' — fl{V') denote the old and new fermionic matrix, A — fl'~^fl and ^ is 
a Gaussian random source vector. The standard deviation of this stochastic estimation can 
be written as 

(7^ = < exp(-2r[A - 1]0 - < exp(-r[A - 1]0 >|*^ 

= det -\2A - 1) - det '^{A). (10) 
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Eq. is valid only if the matrix 2A — 1 is positive definite. If the matrix A has even one 
eigenvalue that is less than or equal to 1/2, the formula in Eq. ^ is not valid, the standard 
deviation is infinite. There is no a priori reason to assume that the matrix A has no small 
eigenvalues. This is a very serious problem that could make the stochastic estimator useless 
in dynamical calculations. If the fermionic matrices Q and fl' are close, i.e. only a few links 
are changed in the update, A = Q'^^Q ^ 1 and consequently det{A) ^ 1 and det(2A — 1) ^ 1 
as well. However for an effective updating method we would like to change the configuration 
at as many links as possible, which makes the occurrence of a small eigenvalue likely. In 
the following we propose a 2-step solution that can always be used to handle the small 
eigenvalues in A. 

First we follow the program of Refs.j^, |^ and replace fl and Q' by reduced matrices 

fi, = l]e-2/(^), = fi'e-2/(^') (11) 

with / a yet to be determined polynomial . We can rewrite eq. ^ as 

^exp(2Tr(/(fi')-/(^^))) 



det det ilr 

= < expi-C[n'-'nr - m >ee exp(2Tr(/((]') - f{n))). (12) 

Since Tr f can be calculated exactly, only the first factor of the last expression is evaluated 
stochastically. Its fiuctuations are minimized if Ar = fl'~^Qr ~ 1- It is difficult to optimize 
the polynomial / both for Q and Q' at the same time, instead we choose / such that ~ 1. 
That also guarantees fi,, ~ 1- Since for staggered fermions the eigenvalues of the matrix 
Q can vary between 4m^ and 16 + 4m^, we choose the polynomial / such the the function 
^2f{x) 1^ is close to one in that range. In practice we use a third order polynomial 

f{x) = ao + Ci2X + a4X^ + a^x^ (13) 

and choose the coefficients ai by minimizing the function 

A = / (-e'^^") - lfpix)dx. (14) 

J4m2 X 

The weight function p{x) should approximate the eigenvalue density distribution of the 
fermionic matrix. We used a linear approximation for p 

p{x) = X, a; G (4m^, 8 + 4m^) 

= 16 + 8m^-x, a; e (8 + 4m^ 16 + 4m^) (15) 



and considered mass values m = 0.01 — 0.1. We have also tried more complex forms for the 
eigenvalue density that included higher order terms, all motivated by free field calculations. 
The results were not very sensitive to the specific choice of p. We do not want to change 
the a parameters depending on the quark mass of the simulation so we decided to use the 
following values in all cases 

ao = -0.34017 
a2 = 0.35645 
a4 = -0.030379 

as = 0.000957. (16) 

The eigenvalues of the reduced matrix fi^ span a smaller range than the original fermionic 
matrix. At mass am = 0.1 the smallest eigenvalue is increased from 4m^ = 0.04 to about 
0.08, while the ratio of the largest to smallest eigenvalue is reduced from (16 + 4m^)/4m^ ^ 
400 to about 14. At am = 0.04 the increase in the smallest eigenvalue is from 0.0064 to 
0.0125, while the reduction in the ratio of the largest to smallest eigenvalue is from 2,500 to 
about 95. 

When expressed in terms of the reduced matrices the acceptance probabilities of Eqs. HJI 
contain a new gauge-action like term 

Pace = min{l, exp(-ASg + 2Af) < expi-CH'^^r - 1]0 >es} 
Pstoch = min{l, exp(-ASg + 2Af) exp{-C[n'-'^r - 1]0} (17) 

with ASg = Sg{V') - Sg{V) aud Af = Tv{f{n') - /(^])). To calculate the expres- 
sion C,*[Q'^^Qr — 1]^ requires multiplications with Q, Q'^^, and with the reduction factors 
exp(— 2/(fi)) and exp(2/(f2')). The exponentials can be expanded in a Taylor series and 
approximated with a few terms. In [0] we found that it was sufficient to keep only 15 terms. 
Later we will argue that it is better to replace Qr and themselves with a finite order 
polynomial. 

To complete the evaluation of the determinant and acceptance probabilities we still have 
to calculate the trace of f{0,). Tr can be expressed as the combination of the plaquette 
and the three 6-link loops of the smeared links V. It is a fairly straightforward calculation 
giving 

Tvf{n) = (-8/34 + 336/?6)EReTrn„-t- 

n 
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12^6[- E Re Trm„ - E Re Tr ^„ + 

n n 

E Re TyJTu - SNtfi E Re Tr Pn] + const, (18) 

n n 

where the summation is over all distinct objects in the lattice. Pn is the Polyakov line that 
gives a contribution on lattices of size = 6. Lattices that are smaller than 5 in any 
direction would have additional contribution of length-6 overlapping loops. Eq. |l^ is not 
valid in that case. The coefficients P are related to the quark mass and the optimized a 



parameters of Eq. [T^ as 

Po = ao + «2 ■ 4m^ + ^4 ■ (4m^)^ + ag ■ (4m^)^ 
P2 = — «2 — 04 ■ 2 ■ 4m^ — ag ■ 3 ■ (4m^)^ 
Pi = 04 + cte ■ 3 ■ 4m^ 

Pe = -ae. (19) 

In Ref. we used a similar reduction using a second order polynomial for /. There 
we determined the 02 and coefficients by trial and error attempting to maximize the 
acceptance rate in Eq. ^. The values we obtained there are consistent with what we would 
determine now with the minimization procedure. Increasing the order of the reduction 
polynomial further gives only slight improvement but would require the evaluation of Tr fl^ 
and higher order terms in Tr /, a considerable computational task. 

The polynomial reduction of the fermionic matrix results in considerable improvement in 
the evaluation of the determinant ratio but it is not sufficient to guarantee that the standard 
deviation of the stochastic estimator is finite or small. To achieve that we now proceed to 
rewrite the reduced fermionic determinant ratio as 

det -^Ar = det -"Aj/" 

= <expi-±^;[Al/--m>^.^^ (20) 
j=i 

with n an arbitrary positive integer. The expectation value is evaluated with n independent 
random vectors and the standard deviation becomes 

= det -"(2^1^" - 1) - det '^(Ar). (21) 

cr^ is finite if none of the eigenvalues of the matrix Ar is smaller than or equal to 2~". This is 
a much easier condition to satisfy then the one before. With the reduced matrix, assuming 



that the smallest eigenvalue of Ar is not smaller than the product of the smallest eigenvalues 
of fir and fl'~^, n = 4 is sufficient to guarantee the finiteness of the standard deviation at 
am = 0.1, and n = 8 is sufficient at am = 0.04. For small standard deviation one might 
need to use larger n values, but with any mass am ^ and matrix A^. it is possible to choose 
n such that the standard deviation is finite and small. The cost of this improvement is that 
we have to evaluate the expression ^*(y4y" — 1)^ n times for one estimate of the determinant. 

The nth root of the matrix A^. can be approximated with polynomials to arbitrary preci- 
sion 1^, |l^]. Since the order of the matrices in the determinant are irrelevant, we write the 
nth root of the matrix as 

= fi' ;l/2n ^1/™ ^ ~V^n_ ^22) 

We found that breaking up Vt' to two identical terms and separating them with the less 
singular Vt^J^ improves the stochastic estimator. The terms Vtl^^ and are approxi- 

mated with polynomials 

Qi/n =fii/nexp(-2/(r])/n) = Qt\fl), (23) 

where P^^^"^ and qI"^ are / and k order polynomials of the fermionic matrices f2 and VL' . To 
reduce the errors of the polynomial approximation we write the exponent in the stochastic 
estimator as 

c[Aii^ - i]e = rp/'"H^^') Qt\^)p!;^''^ mi - cpf'^^\^') Qt'\^') pf'''\n% (24) 

The necessary order for the polynomials P and Q vary with the quark mass but we found 
that in most cases fairly low orders are sufficient. 

At this point the generalization of the partial-global updating method to arbitrary flavor 
number is straightforward. To describe Uf flavors we have to replace the determinant ratio 
in Eq. |^ by its Uf/Ath power. That can be easily done by summing up to only in Eq. 



20|. The polynomials P and Q do not have to be changed and smaller n will be sufficient for 



the same standard deviation of the determinant. 



B. Calculating the determinant ratio 



To illustrate the improved estimator, in this section we calculate the ratio det{Q')/ det(O) 
for a specific configuration pair. We chose {U} from a configuration set that was generated 
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Figure 1: The stochastic estimator for det{Q')/ det{Q) using a) the naive estimator with / = , 



b) the estimator with improved / given in Eqs. |13|, |16| and c) the form of Eq. ^ with n = 8 and 
improved / . 

with the HYP dynamical action at /3 = 5.2, am = 0.1, and SgiV) = 0. The scale at these 
parameter values is ro/a = 3.0(1) and ni-^/mp ^ 0.8 [0]. We updated 300 random links of 
{U} with a heat bath step corresponding to a /3 = 5.2 plaquette gauge action to create the 
{U'} configuration. To calculate the determinant ratio we use Eq. ^ with / = 0, with the 



improved / given in Eqs. |13|, |I6|, and also using the formula of Eq. ^ with n = 8 and the 
improved /. We calculate the expectation value using 500-1500 random vectors. Figure [l| 
shows the stochastic estimator for the three cases. One could not guess from the figure that 
the three estimators describe the same quantity. The naive / = estimator is 15 orders of 
magnitude smaller than the improved ones. How is that possible? The average < e"'^'^^ > 
will be the same for all three estimators, but for the naive one the average will come from 
many almost zero values and an occasional large one. That occasional large value is so rare 
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that we did not even encounter it in 500 samples. The improved estimator with n = 1 looks 
much more reliable and it predicts 

^^0.81(12), n.l. (25) 

The estimator with n = 8 is even better. With only a third of the statistics of the n = 1 
case it predicts the determinant ratio as 

^ = 0.77(2). n = 8. (26) 

That does not mean that the acceptance rate of the partial-global update is close to 80% if 
we update 300 links at a time. The acceptance rate is better described by the expectation 
value 

< min{l, e"^^*} >ee= 0.32(1), n = 1, 

= 0.65(1), n = 8. (27) 

One should remember that the above values correspond to a specific pair of configurations. 
Before calculating the determinant ratio on an ensemble of configurations we first discuss a 
modification of the gauge action. 



C. Choosing the gauge action Sg{V) 

In the previous chapter we showed how to remove the most singular part of the inverse 
fermion matrix by multiplying it with a factor exp(2/(fi)). The change in the fermion 
determinant is compensated by the additional term exp(2 Ti{f {Q') — f (Q))) in the stochastic 
estimator that can be calculated exactly. With f{fl) a third order polynomial Tr f{Q) is 
a combination of the plaquette and 6-link loops of the smeared links as given in Eq. p!8| . 
While it is straightforward to evaluate /(fi), it is not a completely negligible computational 
cost. On the smeared link lattice the plaquette and the 6-link loops are very correlated and 
Tr f{^) can be approximated by the plaquette term only, thus reducing the computational 
overhead. In general we choose the gauge action Sg{V) as 

Sg{V) = 2 Tr fiQ) - 12(3,6n„, E Re Tr P„ - ^ ^ Re Tr Q,. (28) 

Here we have included the Polyakov line term explicitly, i.e. SgiV) is the combination of the 
4 and 6-link gauge loops only, it contains no loops closed because of the periodicity of the 

10 



lattice. Like before, this action should not be used if in any direction the lattice is smaller 



than 5. With this choice for Sg{V) the gauge term in the acceptance probability in Eq. ^ 
simplifies to 

- ASg + 2Af = 12/365^,,6E(ReTrP^ - ReTrP„) + ^ ^(ReTrQ, - ReTrQ,')- (29) 

We can choose the coefficient 7 to account for the 2Tr f{fl) term, or even better, we can 
choose it to maximize the determinant ratios and the acceptance rate. Then we not only 



avoid the computation of the term exp(2 Tr(/(n') — f{^))) in Eq. [l^ but can also increase 
the efficiency of the updating algorithm. 

When we write the gauge action as Sg{U) + SgiV), we break up the gauge term into 
two pieces. We use the first term Sg{U) in the heat bath update and include the second 
one, Sgiy), in the accept-reject term. Such a break-up usually lowers the acceptance rate, 
especially if the second term fiuctuates considerably. With the choice of Eq. ^ Sg{V) 
actually cancels an other fiuctuating term, 2A/, and the algorithm should get more efficient. 
The introduction of the plaquette term proportional to 7 could compromise the efficiency. 
Since the smeared plaquette term does not fiuctuate much, a small 7 coefficient does not 
harm the acceptance rate much. How should we choose the coefficient 7? According to Ref. 
the HYP action with SgiV) = at /3 = 5.2, am = 0.1 has lattice spacing a ~ 0.17fm. In 
the global heat bath update the links of the configurations are updated with a pure gauge 
action of gauge coupling /3 = 5.2. The pure gauge configurations at this coupling are very 
different form the dynamical configurations. The correlation length that characterizes the 
large distance behavior is much smaller on the pure gauge configurations. At short distances 
the average plaquette on the dynamical configurations is < ReTrD >dyn= 1-45 while on the 
pure gauge configurations the average plaquette is much smaller, < ReTrCH > (3=5.2= 1.30. 
The gauge action that we use to create new configurations does not match the dynamical 
action neither at long nor at short distances. One would expect that the partial-global 
update is most effective when the pure gauge configurations of the heat bath step are close 
to the dynamical configurations. This suggests that in order to maximize the efficiency of 
the partial-global update we can try to match the short and/or long distance fluctuations 
of the heat bath and dynamical actions. To match the short distance fluctuations we can 
require that the average plaquette of the heat bath update action and the dynamical action 
are close. This condition requires different 7 coupling at different quark masses and gauge 
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coupling values but choosing 

7 = -0.1 (30) 

offers a good compromise. The choice 7 = 0.0 is not much worse and leads to a somewhat 
simpler action, but in the following we will use 7 = —0.1. By construction now the small 
scale fluctuations of the pure gauge heat bath action and of the dynamical action are about 
the same. The modification also improves the matching of the large distance correlations. 
With the new action, in order to reproduce the lattice spacing a = 0.17fm of the P = 5.2, 
am = 0.1, Sgiy) = action, we have to choose the gauge coupling (3 = 5.65 while keeping 
am = 0.1. The lattice spacing of the pure gauge model at /3 = 5.65 is almost 0.17fm, a very 
good agreement. At smaller quark masses the same lattice spacing requires a smaller gauge 
coupling P suggesting that a slightly larger 7 value would be the optimal one. At this point 
we feel that the difference is not significant to justify a quark mass dependent coupling. 

By choosing the gauge action Sg{V) according to Eq. |28| we modify the dynamical action. 
This modification will not affect the perturbative properties of the smeared link action as the 
terms in Sg{V) are independent of the thin link gauge coupling and will become negligible 
in the continuum limit. At finite lattice spacing the new terms could change the scaling 
behavior of the system and their effect should be investigated. 



IV. DETERMINANT RATIOS WITH THE MODIFIED ACTION 



In this chapter we investigate the fermionic determinant ratios on configurations gener- 



ated with the modified action of Eqs. |2^, |3^. We will use two sets of configurations. Both 
sets contain about 100 8^24 lattices. The first one was created at (3 = 5.65, am = 0.1 and 
has lattice spacing a = 0.17fm (ro/a = 2.95(5)) and m.„/mp ^ 0.70. The second set is 
at couplings (3 = 5.55, am = 0.04 with lattice spacing a = 0.17fm (ro/a = 2.88(6)) and 
m-n/mp ~ 0.55. On both sets we created pairs of configurations by updating a random sub- 
set of the original thin links with a heat bath step corresponding to the thin link pure gauge 
couplings, i.e. [3 = 5.65 and [3 = 5.55. On each pair we calculated the modified determinant 
ratio exp(— AS'g) det~^(fii7'~^) using Eqs. [f|, ^ , ^ with 400(800) random source vectors with 
n = 4(8) break-up of the determinant, i.e. we estimated the determinant value from 100 
independent measurements on each configuration pair. We calculated the determinant both 
with relatively small order polynomials (order 16 to 32) and higher order polynomials (order 
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Figure 2: The distribution of the modified fermionic determinant ratios on configuration set I. a) 
tHB = 3000 hnks are updated with a heat bath step and the determinant ratios are calculated with 
n = 4 (dotted lines) and n = 8 (solid lines) determinant break-up. b) all links of a given direction 
and parity {tHB = 6144) are updated at once and the determinant ratio is calculated with n = 8 
break-up. 

64 to 128) to monitor possible systematical errors. The difference between the small and 
higher order approximations is small and well within the errors of the final results. The 
numbers we present here were obtained with the higher order polynomials. 

Figure shows the distribution of the modified determinant ratios on 80 configuration 
pairs from set I. The histogram of figure ^a corresponds to determinant ratios on con- 
figuration pairs that differ at 3000 links. The solid lines shows the distribution measured 
with n = 8, the dotted lines with n = 4 determinant break-up. The two measurements are 
consistent predicting 



< e"^^« det >seti= 0.80(13), tHB = 3000 



(31) 



for the determinant and 



< min{l, e^"^^^ det -\nn'-')} >seti= 0.54(5), tHB = 3000 



(32) 



for the acceptance rate as defined in Eq. 0. While the n = 4 and n = S measurements 
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Figure 3: The configuration average of the stochastic acceptance rate < -Pstoch >u on the configura- 
tion set I as the function of the number of hnks touched in the heat bath update. Bursts correspond 
to n = 8, octagons to n = 4 break-up of the determinant. The crosses correspond to the acceptance 



agree in their prediction of the determinant ratios, they differ considerably in their standard 
deviation. The average standard deviation as defined in Eq. ^ of the n = 4 calculation is 
a"„=4 = 3.5(8) while for = 8 it is <Jn=8 = 1-8(3). The standard deviation of the determinant 
measurement can influence the autocorrelation time of a simulation as that depends both 
on the effectiveness of the gauge update and on the error of the stochastic estimator. A 
factor of two increase in the standard deviation of the determinant could require up to a 
factor of four increase in the number of stochastic estimators, increasing the autocorrelation 
time accordingly. The extra computational cost of breaking the determinant up to n = 8 
instead of n = 4 parts could be easily compensated with the reduced autocorrelation time. 
Whether it is worth using even larger number of terms should be investigated at different 
quark masses separately. 

With the heat bath update we change a random set of links of given direction and parity. 
On an 8^24 conflguration a maximum of 6144 links can be changed at once. In flgure H/b 
we show the modified determinant ratio distribution when we update all 6144 links of a 



rate using the determinants as defined in Eq. ^ 
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randomly chosen direction and parity. This result was obtained with n = 8 break-up. The 
average of the determinant ratios 

< e-^^« det -\nn'-^) >soti= 0.71(10), tnB = 6144, (33) 

and the acceptance rate as defined in Eq. |^ 

< min{l, e"^^^ det -\nn'-^)} >seti= 0.44(4), tnB = 6144 (34) 

are not much diff'erent from the previous tnB = 3000 values, though the average standard 
deviation is worse, cr„=8 = 2.7(6). To us this is a surprising result: on a lOfm^, 8^24 lattice we 
can perform a heat bath update on all the links in a given direction and parity, the maximum 
that can be updated on this volume simultaneously, and accept this change with close to 
50% probability. The configuration average of the stochastic acceptance rate < Pgtoch >u 
of Eq. H is not that high. Figure || compares the average stochastic acceptance rate as the 
function of the links touched in the heat bath update both for n = 4 and n = 8 determinant 
break-up and the acceptance rate from the determinant as defined in Eq. |. With n = 8 
the stochastic acceptance rate is close to 20% if ins = 6144 and about 30% if ins = 3000. 
The stochastic acceptance rate with n = 4 is somewhat lower. Even though these values 
are smaller than the maximal ones predicted by the determinants themselves, they are still 
quite large. What parameters would provide the best choice in an actual simulation depends 
on many things: the number of links that effectively change in an update step, the cost of 
increasing the breakup of the determinant, and on the autocorrelation time of the simulation. 
The study of these questions is beyond the scope of the present paper and we will return to 
them in a future publication. 

Figures || and || are the same as figures ||/a and || but for configuration set II. Most 
fermionic methods lose some efficiency at smaller quark masses and the stochastic estimator 
is no exception. While the average of the determinant ratios 

< e"^^" det -\nn'-^) >setii= 0.8(2), tnB = 3000, (35) 

is not much different from the set I configurations, the acceptance rate as defined in Eq. || 
is smaller 

< min{l,e-'^^^det -^(^]^]'-^)} >sctii= 0.35(7), tne = 3000. (36) 



15 



50 



: tHB=3000 

40 =1 



30 



20 - 



10 - 




2 4 6 

e-^^«det-^(A) 

Figure 4: The distribution of the modified fermionic determinant ratio on configuration set II. 
tHB = 3000 hnks are updated with a heat bath step and the determinant ratios are calculated with 
n = 4 (dotted lines) and n = 8 (solid lines) break-up. 

Whether this decrease is due to the smaller quark mass or reflects the fact that the pure gauge 
heat bath action does not match the dynamical action well is worth further investigation. 
To match the stochastic acceptance rate of set I with n = A determinant break-up we have 
to use n = 8 on set II as flgure |5| shows. 

V. CONCLUSION 

In this paper we proposed an improved method to calculate the fermionic determinant 
of dynamical conflgurations. The method is very general but relies on the smoothness of 
smeared gauge links. To test the method we considered dynamical conflgurations, updated a 
large subset of their links with a pure gauge heat bath step, and calculated the ratios of the 
fermionic determinants on the old and new conflgurations. We found that even if all the links 
of a given direction and parity of an 8^24, lOfm^, mp/m-,, = 0.7 conflguration are updated 
at once, (6144 links in all), the fermionic determinant ratio is still fairly large and such a 
change would be accepted by a Metropolis accept-reject step with about 50% probability. 
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Figure 5: The configuration average of the stochastic acceptance rate < -Pstoch >u on the con- 
figuration set II as the function of the number of hnks touched in the heat bath update. Bursts 
correspond to n = 8, octagons to n = 4 break-up of the determinant. The cross corresponds to the 
acceptance rate using the determinants as defined in Eq. ^. 

Using only a single stochastic estimator for the determinant reduces the acceptance rate 
to 20% but still offers an effective update. On configurations with smaller quark masses 
the stochastic estimator loses some efficiency. When mp/m^ = 0.55 only about half that 
many links can be updated at one time vi^ith 20% stochastic acceptance rate though the 
determinants stay about the same as v^ith larger quark masses. 

We have not used the fully improved method in dynamical simulations yet, nor did we 
optimize all its parameters. The optimization requires tuning the parameters of the action 
and calculating autocorrelation times with different determinant break-up and updating 
steps. This work is in progress and the results will be reported in a forthcoming publication. 
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